Sampling Challenges
According to Andrew Gelman (2008)
When you have computational problems, often there’s a problem with your model.
When sampling goes poorly, often there is a problem with the model:
See lecture 2-4.
[1] 44.88692 48.25894 48.44770 50.16844 49.55441 49.64544 48.49286 48.91631
[9] 50.18508 53.58071 51.02436 50.97705 50.22559 52.74691 51.13053 52.56043
[17] 52.78518 51.04284 51.18040 50.38051 49.88772 49.00229 47.69577 52.15505
[25] 54.26835 49.64690 49.63739 52.25737 47.74531 51.14175 52.37880 50.37451
[33] 50.81945 52.62929 49.75200 50.66841 52.09346 45.51304 50.30989 48.95920
[41] 49.78614 48.94043 49.87016 52.75281 48.30467 49.73333 50.92056 52.59823
[49] 52.43258 50.08343 47.82025 46.46541 55.68929 50.39334 49.46525 49.26600
[57] 48.82329 49.49899 46.46622 48.76405 50.32061 49.12468 52.92599 51.42727
[65] 51.01789 53.63535 49.35317 52.96504 50.06350 49.25569 49.72310 49.79884
[73] 49.89066 49.14788 53.12937 49.11339 48.83087 47.81584 47.54861 46.55724
[1] 50.16057
[1] 2.002828
brm() modelbrm summary Family: gaussian
Links: mu = identity; sigma = identity
Formula: y ~ 1
Data: list(y = y) (Number of observations: 80)
Draws: 4 chains, each with iter = 10000; warmup = 5000; thin = 1;
total post-warmup draws = 20000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 50.16 0.23 49.71 50.60 1.00 17282 13194
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 2.03 0.16 1.74 2.38 1.00 17442 13195
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
\[ Mass = \frac{a_1}{1 + \exp (-b_1 (age - c_1))} + \frac{a_2}{1 + \exp (-b_2 (age - c_2))}\]
brm modelpriors <-
prior(normal(9, 5), nlpar = "a1") +
prior(normal(4, 5), nlpar = "a2") +
prior(normal(0.5, 0.001), nlpar = "b1") +
prior(normal(0.1, 0.001), nlpar = "b2") +
prior(normal(10, 5), nlpar = "c1") +
prior(normal(17, 5), nlpar = "c2")
fm <- brm(
bf(Mass ~ a1 / (1 + exp(-1 * b1 * (Age - c1))) +
a2 / (1 + exp(-1 * b2 * (Age - c2))),
a1 + b1 + c1 + a2 + b2 + c2 ~ 1,
nl = TRUE),
data = D,
prior = priors,
backend = "cmdstan",
chains = 4,
cores = 4,
iter = 1e4,
refresh = 2000,
seed = 347922,
save_pars = save_pars(all = TRUE)
)Warning: 2070 of 20000 (10.0%) transitions ended with a divergence.
See https://mc-stan.org/misc/warnings for details.
brm summary Family: gaussian
Links: mu = identity; sigma = identity
Formula: Mass ~ a1/(1 + exp(-1 * b1 * (Age - c1))) + a2/(1 + exp(-1 * b2 * (Age - c2)))
a1 ~ 1
b1 ~ 1
c1 ~ 1
a2 ~ 1
b2 ~ 1
c2 ~ 1
Data: D (Number of observations: 6)
Draws: 4 chains, each with iter = 10000; warmup = 5000; thin = 1;
total post-warmup draws = 20000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
a1_Intercept 6.70 2.02 2.91 11.00 1.02 864 1227
b1_Intercept 0.50 0.00 0.50 0.50 1.00 1322 1644
c1_Intercept 14.62 1.55 11.13 16.86 1.01 1528 2413
a2_Intercept 7.52 3.41 0.18 13.72 1.01 1054 1192
b2_Intercept 0.10 0.00 0.10 0.10 1.00 2413 3242
c2_Intercept 19.15 4.70 9.63 28.04 1.00 1112 2808
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.55 0.38 0.19 1.58 1.01 259 118
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
brm modelpriors <-
prior(normal(9, 5), nlpar = "a1") +
prior(normal(4, 5), nlpar = "a2") +
prior(normal(0, 2), nlpar = "b1") +
prior(normal(0, 2), nlpar = "b2") +
prior(normal(10, 5), nlpar = "c1") +
prior(normal(17, 5), nlpar = "c2")
fm <- brm(
bf(Mass ~ a1 / (1 + exp(-1 * b1 * (Age - c1))) +
a2 / (1 + exp(-1 * b2 * (Age - c2))),
a1 + b1 + c1 + a2 + b2 + c2 ~ 1,
nl = TRUE),
data = D,
prior = priors,
backend = "cmdstan",
chains = 4,
cores = 4,
iter = 1e4,
refresh = 2000,
seed = 984575,
save_pars = save_pars(all = TRUE)
)brm summary Family: gaussian
Links: mu = identity; sigma = identity
Formula: Mass ~ a1/(1 + exp(-1 * b1 * (Age - c1))) + a2/(1 + exp(-1 * b2 * (Age - c2)))
a1 ~ 1
b1 ~ 1
c1 ~ 1
a2 ~ 1
b2 ~ 1
c2 ~ 1
Data: D (Number of observations: 6)
Draws: 4 chains, each with iter = 10000; warmup = 5000; thin = 1;
total post-warmup draws = 20000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
a1_Intercept 9.73 3.75 2.48 16.65 1.18 17 30
b1_Intercept 0.28 1.51 -4.82 2.81 1.14 38 35
c1_Intercept 11.93 3.93 2.68 17.83 1.13 34 102
a2_Intercept 4.56 5.36 -5.17 13.48 1.24 12 125
b2_Intercept 0.44 1.27 -2.69 3.17 1.09 69 233
c2_Intercept 17.41 4.33 10.22 26.17 1.09 33 525
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.60 0.46 0.09 1.85 1.07 53 31
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
Warning messages:
1: Parts of the model have not converged (some Rhats are > 1.05). Be careful when analysing the results! We recommend running more iterations and/or setting stronger priors.
2: There were 3227 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
When stan models fail, they often fail badly and with a lot of warning messages.
control = list(adapt_delta = 0.99,
max_treedepth = 15)